join.contigs <- function(contig,blast.data,contig.info,out.file=NULL){
	#tries to pick series of contigs corresponding to blast hit of 
	#EST versus 454 contigs
	#
	#contig is the focal contig
	#blast.data is the output from 454 vs ests
	#contig.info is length and no of reads per contig
	require(GeneR)
	
	tmp.blast <- blast.cDNA2[blast.data$Query.id %in% contig,]
	if(dim(tmp.blast)[1]==0) {
		cat('no blast hit found for ',contig,'\n')
		return(0)
		}
	tmp.blast <- tmp.blast[order(tmp.blast$bit.score,decreasing=TRUE),]
	#think about adding another column for number of reads

	tmp.inc <- rep(FALSE,dim(tmp.blast)[1])
	max.q <- max(c(tmp.blast$q.start,tmp.blast$q.end))
	tmp.olap <- rep(0,max.q)

	#picks best contigs, starting from the best bit score and avoiding big overlaps
	for(i in 1:dim(tmp.blast)[1]){
		if(sum(tmp.olap[tmp.blast[i,"q.start"]:tmp.blast[i,"q.end"]])>8) next
		tmp.olap[tmp.blast[i,"q.start"]:tmp.blast[i,"q.end"]] <- 1
		tmp.inc[i] <- TRUE
		}
	
	tmp.blast2 <- tmp.blast[tmp.inc,]
	tmp.sc.contigs <- unique(tmp.blast2$Subject.id[order(tmp.blast2$q.start)])

	if(is.null(out.file)){
		tmp.fileout <- paste("test_files/",tmp.blast2$Query.id[1],"_matched.embl",sep="")
		}else{
		tmp.fileout <- out.file
		}
	
	writeEmblLine(file=tmp.fileout,code="ID",append=F,text=paste(tmp.blast2$Query.id[1],"_matched",sep=""))

	tmp.sc.contigs2 <- sub('lcl\\|','',tmp.sc.contigs)

	tmp.pst <- paste(tmp.sc.contigs2,collapse="; ")
	tmp.cont.info <- sc_contig_info2[match(tmp.sc.contigs2,sc_contig_info2$sc.contig),]
	tmp.cont.info$length <- as.numeric(tmp.cont.info$length)
	tmp.cont.info$numreads <- as.numeric(tmp.cont.info$numreads)
	tmp.cont.info$direction <- factor("forward",levels=c("forward","reverse"))
	writeEmblLine(file=tmp.fileout,code="CC",text=tmp.pst)
	writeEmblLine(file=tmp.fileout,code="FT",header="source", 
		text=paste("1",sum(tmp.cont.info$length)+100*(dim(tmp.cont.info)[1]-1),sep=".."),nextfield=F)
	writeEmblLine(file=tmp.fileout,code="FT",
		text="/organism=\"Lagopus lagopus\"",nextfield=F)
	writeEmblLine(file=tmp.fileout,code="FT",
		text="/mol_type=\"genomic DNA\"",nextfield=F)

	#say where the contigs are
	tmp.l <- 1
	tmp.cont.info$concat.start <- 0
	tmp.cont.info$concat.end <- 0

	#begin big loop to annotate each contig

	for(i in 1:length(tmp.sc.contigs)){
		tmp.cont.info$concat.start[i] <- tmp.l
		tmp.l <- tmp.l + tmp.cont.info$length[i] +100 -1
		tmp.cont.info$concat.end[i] <- tmp.l
		tmp.l <- tmp.l +1
		writeEmblLine(file=tmp.fileout,code="FT",header="misc_feature",
			text=paste(tmp.cont.info$concat.start[i],"..",tmp.cont.info$concat.end[i]-100,sep=""),nextfield=F)
		#change header to misc_feature?
		writeEmblLine(file=tmp.fileout,code="FT",text=paste("/note=\"",tmp.cont.info$sc.contig[i],"\"",sep=""),nextfield=F)
		tmp.blast.cont <- tmp.blast2[tmp.blast2$Subject.id==tmp.sc.contigs[i],]
		tmp.seqi <- strReadFasta("grouse_seq_cap/AllGrouseContigs.fasta",name=tmp.cont.info$sc.contig[i])
	
		if(tmp.blast.cont$s.start[1]>tmp.blast.cont$s.end[1]){ #reverse complement
			tmp.seqi <- strComp(tmp.seqi)
			writeEmblLine(file=tmp.fileout,code="FT",text="/note=\"reverse\"",nextfield=F)
			tmp.cont.info$direction[i] <- "reverse"
			}else{
			writeEmblLine(file=tmp.fileout,code="FT",text="/note=\"forward\"",nextfield=F)
			}
		if(i < length(tmp.sc.contigs)){
			writeEmblLine(file=tmp.fileout,code="FT",header="gap",
			text=paste(tmp.cont.info$concat.end[i]-100+1,"..",tmp.cont.info$concat.end[i]-1,sep=""),nextfield=F)
		writeEmblLine(file=tmp.fileout,code="FT",header="",
			text="/estimated_length=unknown",nextfield=F)
			}
		
		#add markup for blast hits
		if(tmp.blast.cont$s.start[1]>tmp.blast.cont$s.end[1]) {
			tmp.s.start <- tmp.cont.info$length[i]+1 - tmp.blast.cont$s.end
			tmp.s.end <- tmp.cont.info$length[i] - tmp.blast.cont$s.start +1
			}else{
			tmp.s.end <- tmp.blast.cont$s.end
			tmp.s.start <- tmp.blast.cont$s.start
			}
		tmp.s.end2 <- tmp.cont.info$concat.start[i] + tmp.s.end-1
		tmp.s.start2 <- tmp.cont.info$concat.start[i] + tmp.s.start-1
	
		tmp.tf <- FALSE
	
		for(j in 1:length(tmp.s.start2)){
			if(i == length(tmp.sc.contigs)&j == length(tmp.s.start2)) tmp.tf <- TRUE
			writeEmblLine(file=tmp.fileout,code="FT",header="misc_feature",
				text=paste(tmp.s.start2[j],"..",tmp.s.end2[j],sep=""),nextfield=F)
			#change header to misc_feature?
		
			writeEmblLine(file=tmp.fileout,code="FT",
				text=paste("/note=\"blast hit to query: q.start, ",tmp.blast.cont$q.start[j], 
				"; q.end, ",tmp.blast.cont$q.end[j],
				"; alignment length, ",tmp.blast.cont$alignment.length[j],
				"; mismatches, ",tmp.blast.cont$mismatches[j],
				"; gap openings, ",tmp.blast.cont$gap.openings[j],
				"; bit.score, ",tmp.blast.cont$bit.score[j],"\"",sep=""), nextfield=tmp.tf)
			}
	
		if(i==1){
			tmp.seqout <- tmp.seqi
			}else{
			tmp.seqout <- paste(tmp.seqout,paste(rep("N",100),collapse=""),tmp.seqi,sep="")
			}
	
	
		}
	tmp.fileoutf <- sub('\\.[a-z]*','\\.fasta',tmp.fileout)
	placeString(tmp.seqout,seqno=0)
	setStrand(0)
	writeEmblSeq(tmp.fileout)
	
	if(tmp.fileoutf!=tmp.fileout){
		writeFasta(file=tmp.fileoutf,name=sub('\\.[a-z]','',tmp.fileout))
		cat("Sequence written to ",tmp.fileout," and ",tmp.fileoutf,"\n")
		}else{
		cat("Sequence written to ",tmp.fileout,"\n")
		}
	placeString("",seqno=0)
	tmp.cont.info
	}
#join.contigs("cDNA_1405-1",blast.data=blast.cDNA2,contig.info=sc_contig_info2,out.file="test_files/c1405tst.embl")